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Detailed protocol 


#Usage: pythons .py strain_list 
#strain_list format: ERR137225_1.fastq.gz --> ERR137225 
import sys 
f1=open(sys.argy[1]).readlines() 
for strain in f1: 
strain=strain.rstrip() 
step1="/home/zty/sickle-master/sickle pe -t sanger -f"+straint"_1.fastq.gz -r"+straint"_2.fastq.gz-o "+straint+"_1.fastq -p "+strain+"_2 fastq -s 
"+straint"_s.fastq"+"\n" 
step2="/home/zty/bowtie2-2.2.9/bowtie2 -x /home/zty/tb_h37rv.fasta -1 "+straint"_1 .fastq -2 "+straint"_2.fastq -U "+strain+"_s.fastq -S "+straint".sam"+"\n" 
step3="samtools view -bhSt /home/zty/tb_h37rv.fasta.fai "+straint".sam -o "+straint+".paired.bam"+"\n" 
step4="samtools sort "+straint+".paired.bam -o "+straint".sort.bam"+"\n" 
step5=("depth=$( samtools depth "+ strain+".sort.bam | awk '{s+=$3}END {print s/NR}' )"+"\n" 
"coverage=$( samtools depth "+strain+".sort.bam | awk 'END{print NR/4411532}' )"+"\n" 
"a=$(($( echo $depth | awk ‘{printf (\"%.f\",$1)}' )))"+"\n" 
"if [\"$a\" -ge 10 ] && (echo ${coverage} 0.95 | awk '!($1>=$2){exit 1}');then"+"\n" 
"\t't"samtools mpileup -q 30 -Q 30 -Bf /home/zty/tb_h37rv.fasta "+straint+".sort.bam > "+straint".pileup"+"\n" 
"\t'+"b=$(($( echo $depth | awk ‘{printf (\"%.f\" $1 )}' /10))"+"\n" 
"\t'+"if $b < 5;then"+"\n" 
"\t'+"\t"+"java -jar /home/zty/VarScan.v2.3.9.jar mpileup2snp "+strain+".pileup --min-coverage 5 --min-reads2 2 --min-avg-qual 30 --min-var-freq 0.75 --p-value 
99e-02 > "+strain+".varscan"+"\n" 
"\t't+"else"+"\n" 
"\t'+"\t"+"java -jar /nome/zty/VarScan.v2.3.9 jar mpileup2snp "+straint" pileup --min-coverage $b --min-reads2 2 --min-avg-qual 30 --min-var-freq 0.75 --p- 
value 99e-02 > "+straint".varscan"+"\n" 
"\t +" fi" +"\n" 
"\t't"java -jar /home/zty/VarScan.v2.3.9.jar mpileup2cns "+strain+" pileup --min-coverage 3 --min-avg-qual 20 --min-var-freq 0.75 --strand-filter 0 --min-reads2 
2 >"+straint".cns"+"\n" 
"\t'+"awk -F 'L]' '{if($9==0 || $10==0)$0=\"\"selse print $0} "+strain+".varscan > "+straint".vars"+"\n" 
"\t'+"perl /home/zty/script/varscan_work_flow/0.1_PE_IS_filt_Rv.pl /nome/zty/script/varscan_work_flow/PPE_INS_list/PPE_INS_loci_Rv.list "+strain+".vars > 
"+straint".var.ppe"+"\n" 


"\t'"+"perl /home/zty/script/varscan_work_flow/1_format_trans.pl"+straint+".var.ppe > "+strain+".var.for"+"\n" 
"\t"+"cut -f2,3,4 "+straint+".var.for > "+straint".snp"+"\n" 


"\t'+"rm "+straint".sam "+straint+".varscan "+straint".paired.bam "+strain+"_s.fastq "+strain+"_1.fastq "+straint"_2.fastq "+strain+".var.for "+straint+".var.ppe 
"+straint" pileup"+"\n" 

"else"+"\n" 

"\t'+"echo "+straint" do not meet criteria: ${depth} "+"${coverage} >> discard"+"\n" 


"\t'+"rm "+straint+".sam "+straint+".varscan "+straint".paired.bam "+strain+"_s.fastq "+strain+"_1.fastq "+straint"_2.fastq "+strain+".var.for "+straint".var.ppe 
"+straint" pileup"+"\n" 

"fi"+"\n" 

) 


open("pair_end.sh","a").write(step1+step2+step3+step4+step5) 
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